Evaluating the Performance of Single and Double Moment Microphysics Schemes 
During a Synoptic-Scale Snowfall Event 

Andrew L. Molthan 

NASA Short-term Prediction Research and Transition (SPoRT) Center 
NASA Marshall Space Flight Center, Huntsville, AL 

24 th Conference on Weather and Forecasting / 20 th Conference on Numerical Weather Prediction 

Increases in computing resources have allowed for the utilization of high-resolution weather forecast 
models capable of resolving cloud microphysical and precipitation processes among varying numbers of 
hydrometeor categories. Several microphysics schemes are currently available within the Weather 
Research and Forecasting (WRF) model, ranging from single-moment predictions of precipitation 
content to double-moment predictions that include a prediction of particle number concentrations. 

Each scheme incorporates several assumptions related to the size distribution, shape, and fall speed 
relationships of ice crystals in order to simulate cold-cloud processes and resulting precipitation. Field 
campaign data offer a means of evaluating the assumptions present within each scheme. 

The Canadian CloudSat/CALIPSO Validation Project (C3VP) represented collaboration among the 
CloudSat, CALIPSO, and NASA Global Precipitation Measurement mission communities, to observe cold 
season precipitation processes relevant to forecast model evaluation and the eventual development of 
satellite retrievals of cloud properties and precipitation rates. During the C3VP campaign, widespread 
snowfall occurred on 22 January 2007, sampled by aircraft and surface instrumentation that provided 
particle size distributions, ice water content, and fall speed estimations along with traditional surface 
measurements of temperature and precipitation. In this study, four single-moment and two double- 
moment microphysics schemes were utilized to generate hypothetical WRF forecasts of the event, with 
C3VP data used in evaluation of their varying assumptions. Schemes that incorporate flexibility in size 
distribution parameters and density assumptions are shown to be preferable to fixed constants, and 
that a double-moment representation of the snow category may be beneficial when representing the 
effects of aggregation. These results may guide forecast centers in optimal configurations of their 
forecast models for winter weather and identify best practices present within these various schemes. 
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1. Introduction 

Increases in computing power have lead to the experimen- 
tal or operational use of high resolution weather forecast 
models that attempt to explicitly resolve precipitation pro- 
cesses by using bulk-water microphysics schemes capable 
of predicting both the mass content and size distribution of 
several hydrometeor species. These approaches have been 
used to predict the convective mode for severe weather events 
(Kain et al. 2006) and the development of mesoscale snow 
bands responsible for heavy snowfall (Bernardet et al. 2008). 

Several bulk water microphysics schemes are available 
within the Weather Research and Forecasting (WRF) model 
as of version 3.1.1, with varying numbers of simulated 
hydrometeor classes and methods for estimating their size 
distributions, densities, and fall speeds. In order to evaluate 
these various assumptions, field campaign data are necessary, 
providing measurement of particle size distributions, ice or 
liquid water content, inference of particle bulk densities, and 
their terminal fall speeds. Toward this goal, the Canadian 
CloudSat/CALIPSO Validation Project (C3VP) sought to 
obtain in situ observations of ice crystals and aggregates 
in order to evaluate various microphysics schemes, and to 
serve as a basis for evaluating satellite retrievals of cloud 
properties from current sensors and future members of the 
NASA Global Precipitation Measurement (GPM) mission. 
Centered in southern Ontario and managed at the Canadian 
Centre for Atmospheric Research Experiments (CARE), the 
C3VP campaign provided aircraft and surface measurements 
of ice crystals, dual-polarimetric radar data, and traditional 
surface weather observations during a synoptic- scale snow- 
fall event on 22 January 2007. Herein, discussion focuses on 
the evaluation of WRF model forecasts for the 22 January 
2007 event, utilizing aircraft and surface observations to 
evaluate the assumptions and overall performance of several 
bulk water microphysics schemes currently available to the 
operational forecasting community. 

2. Characteristics of Single and Double 
Moment Bulk Water Microphysics Schemes 

Single- and double-moment schemes vary by particle 
size distribution function, methods of calculating distribution 
shape parameters, relationships between mass and diameter, 
relationships between diameter and terminal fall speed, and 
a variety of assumptions within simulated microphysical 
processes. With the exception of the Thompson scheme, all 
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single- or double-moment microphysics scheme used in this 
study use a form of the gamma distribution: 

N X (D) = N 0X D^e-^\ (1) 

where N ox is referred to as the size distribution intercept, 
p x is the dispersion parameter, and A x is the slope parameter. 
In the following analyses, the subscript v is replaced with 
4 s’ to denote references to the snow category. Marshall and 
Palmer (1948) determined that populations of large, precip- 
itating ice crystals could be represented as an exponential 
size distribution, a special case of the gamma distribution 
(1) where p s is set to zero: 

N s (D)=N os e~^ D (2) 

The total mass content within the size distribution can 
be determined by integrating the product of (1) or (2) and 
a mass-diameter relationship. Locatelli and Hobbs (1974) 
determined mass-diameter, or M(D) relationships for several 
types of crystals in a power-law form: 

M{D) = a m D bm (3) 

Locatelli and Hobbs (1974) observed that the coefficient 
a m varied by crystal type and degree of riming, while the 
exponent values of b m averaged near 2 and suggested that 
the mass of a crystal was proportional to cross-sectional area. 
If crystals are represented by an equivalent diameter sphere 
with an effective density of p s , then the coefficient a m = 
|p 5 and exponent b m = 3. By combining (2) and (3), and 
integrating over the entire size distribution the total mass 
content can be determined: 


M= f a w D bm N n ,D^e- XsD dD = 
Jo 


a m N os T(l H - /L T~ ^m) 

n l + H s +b m 


( 4 ) 

Bulk water microphysics schemes are categorized in terms 
of the number of predicted moments, or M n , of the size 
distribution. The moment of a size distribution is a statistical 
property, the integrated product of the diameter raised to the 
power n and the number concentration of the same diameter. 
In terms of the gamma size distribution, the n th moment is 
defined as: 


M n = 


N 0 sT( 1 + p s + n) 

24 +gs+n 


( 5 ) 


Following this terminology, a single-moment microphysics 
scheme predicts one moment, the mass content (or M^ m ) of 
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the particle size distribution for each precipitating species. 
By predicting the total mass through simulated microphysics 
processes, remaining terms can be determined by assigning 
fixed values or functions to restrict remaining parameters. As 
an example, an assignment of N os and M(D) allow for the 
calculation of X s based upon the predicted snow mass content 
(paQs) acquired from simulated microphysical processes: 

fa — f + Ps + b m ) \ l +Vs+bm 

\ PaQs ) 

3. The 22 January 2007 Snowfall Event and 
Available C3VP Data Sets 

Moderate to heavy snowfall occurred over southern On- 
tario in advance of a warm frontal boundary, associated with 
a midlatitude cyclone that traveled along the U.S. -Canadian 
border on 22 January 2007. Precipitation began at the CARE 
site around 0200 UTC and continued through 0800 UTC, 
with the bulk of the precipitation occurring between 0600 and 
0800 UTC, and a liquid equivalent total of 2.8 mm. Surface 
temperatures during the same period hovered near -9°C. The 
broad shield of warm frontal cloud cover and precipitation 
was sampled by the CloudSat radar, an instrumented Convair- 
580 aircraft, and the operational, C-band, dual-polarimetric 
radar at King City, Ontario. The Convair-580 was equipped to 
measure temperature, relative humidity, hydrometeor content 
with a counterflow virtual impactor (CVI, Twohy et al. 1997), 
and particle size distributions (PSDs) via Particle Measuring 
Systems (PMS) 2D-P and 2D-C probes. 

In order to evaluate forecast model performance, use of 
aircraft data herein focuses on two portions of the flight track 
that represent complete vertical profiles: a descending, non- 
Lagrangian spiral obtained near the site of the King City 
radar, and the ascending departure on a southeast heading 
(Fig. 1). Imagery from the 2D-P and 2D-C probes were used 
to construct PSDs at five second increments of flight time, 
provided by A. Heymsfield of NCAR, with methods applied 
to avoid the adverse effects of small particles resulting 
from the shattering of large crystals on the probe housings 
(Heymsfield et al. 2008). The first and second moments were 
used to estimate the intercept and slope of exponential size 
distributions fit to each PSD (Heymsfield et al. 2004). Each 
PSD is accompanied by a measurement of ice water content 
provided by the CVI. By distributing the CVI estimate 
of total ice mass among the equivalent diameter spheres 
within the PSD, an estimate of the effective bulk density 
(Heymsfield et al. 2004) can be obtained and compared 
against forecast model assumptions. 

4. Generation of WRF Model Forecasts 

Comparisons between model performance and field cam- 
paign measurements require a plausible forecast of the event. 
Shi et al. (2010) reproduced the characteristics of the event 
using a triply nested, 9-3-1 km WRF model domain con- 
figuration, which was subsequently used by Molthan et al. 
(2010) to evaluate the assumptions of the NASA Goddard 
six-class, single-moment, bulk water microphysics scheme. 



Fig. 1. Overview of some observational datasets used herein and available 
during the C3VP campaign. Aircraft profiles used in this study are color- 
coded to represent the descending spiral (red) and ascending departure 
(blue), and repeated in subsequent figures. The crosshair represents the 
location of the dual-polarimetric, C-band radar at King City, Ontario, with 
range rings at 50 km intervals. The dashed line to the northwest represents 
radar cross-sections obtained at the 331° azimuth in the direction of the 
CARE site. 


Here, the configuration of Molthan et al. (2010) (Table 1) 
was used to generate several additional forecasts of the 
22 January 2007 event, modifying the original configuration 
to permit experiments using a variety of single- and double- 
moment schemes available within the Advanced Research 
(ARW) version of the WRF model as of version 3.1.1, and 
the opportunity to evaluate a new scheme proposed by Fin 
and Colle (2010). In total, six forecasts were generated to 
evaluate each scheme’s ability to reproduce aircraft measure- 
ments of temperature, relative humidity, and properties of ice 
crystal size distributions, in addition to radar reflectivity and 
liquid equivalent precipitation. 

a. Microphysics Schemes Used to Simulate the Event 

The schemes used in this study offer a variety of methods 
to determine values of A*, N os , mass-diameter, and diameter- 
fall speed relationships. Six schemes were investigated, in- 
cluding both single and double moment predictions, and 
some characteristics of each are described here. The reader 
is strongly encouraged to review the cited references for 
additional information beyond the size distribution param- 
eterizations described here. The Goddard scheme adopts 
the methodology of Fin et al. (1983), assigning a fixed 
value for N os and a spherical shape representation where 
the effective bulk density of snow crystal populations is 
fixed. The WRF six-class, single-moment (WSM6) and the 
WRF six-class, double-moment (WDM6) schemes assume 
a spherical shape for snow crystals and a fixed, effective 
bulk density is used to define M(D), but the distribution 
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Table 1. Configuration of the WRF model for simulation of the 22 January 2007 snowfall event, following Shi et al. (2010) and Molthan et al. (2010). 


Physical Process 

Parameterization Scheme 

Notes 

Boundary Layer 

Mellor-Yamada-Janjic 

Janjic (1990, 1996, 2002) 

Longwave Radiation 

Rapid Radiative Transfer 

Mlawer et al. (1997) 

Shortwave Radiation 

Dudhia Scheme 

Dudhia (1989) 

Land Surface Processes 

NOAH Land Surface Model 

Ek et al. (2003) 

9 km Cumulus Parameterization 

Grell-Devenyi Scheme 

Grell and Devenyi (2002) 

3,1 km Microphysics Parameterization 

Varied 

See text 


intercept N os is determined as a function of temperature 
following Houze et al. (1979). The WDM6 is equivalent to 
the WSM6 except that the rain category is predicted with 
two moments, both the total mass and number concentration. 
The Stony Brook University scheme by Lin and Colle 
(2010) (SBU-Lin) uses a temperature dependent relationship 
for N os (T ) based upon Houze et al. (1979), along with a 
mass-diameter and diameter-fall speed relationship that is 
determined from a diagnosed riming factor, or Rj (Lin et al. 
2010). The Thompson scheme differs from other single- 
moment schemes by assuming size distribution combin- 
ing exponential and gamma shapes, a non-spherical mass- 
diameter relationship, and predicting additional moments 
of the size distribution based upon temperature-dependent 
relationships between M 2 and other M n acquired from aircraft 
field campaign measurements (Field and Heymsfield 2003; 
Thompson et al. 2008). 

Double-moment schemes predict an additional moment, 
and provide additional information to better define the size 
distribution parameters of N os and X s . In the Morrison 
scheme, the only scheme evaluated in this study that in- 
cludes a double-moment representation of ice, both the 
mass and number concentration are predicted, based upon 
gamma size distributions for each hydrometeor class. Among 
these schemes, the Goddard, WSM6, WDM6 and Thompson 
schemes include prediction of the graupel class. The SBU- 
Lin scheme does not include a separate graupel class but 
instead incorporates variable characteristics of snow crystals 
dependent upon their degree of riming. The Morrison scheme 
does not include the prediction of graupel. Prediction of grau- 
pel is a substantial difference among the schemes presented 
here since the additional class would provide a variety of 
sources and sinks related to the production of an entirely 
separate category, however, observations for the 22 Jan- 
uary 2007 event suggest that snow crystals and aggregates 
were the overwhelming particle type (Petersen et al. 2007), 
and that the simulation of graupel is not key to reproducing 
the character of the event. Selected characteristics of each 
scheme and their relevant parameters are listed in Table 2 
and Table 3, respectively. 

5. Surface Temperature and Precipitation 

Throughout the forecast period, each of the single- or 
double-moment microphysics schemes produced a unique 
depiction of storm total precipitation, related to their intrinsic 
assumptions and simulated processes (Fig. 2). Results herein 
focus on comparisons between model outputs and C3VP 
campaign data to determine which of the available schemes 


would be most appropriate for this specific event. Broader 
assessments of the relative strengths and weaknesses of each 
scheme would require a simliar analysis over multiple events, 
and a validation campaign encompassing other geographic 
regions and types of events. 

Comparisons between CARE site air temperatures and 
the two meter temperatures from the nearest grid point of 
each forecast are shown in Fig. 3a, along with storm total, 
liquid equivalent accumulations of precipitation. Modeled 
surface temperatures were initialized with an apparent warm 
bias obtained from GFS initial conditions, transitioned to 
a slight cool bias beginning around 1800 UTC on 21 Jan- 
uary, and then closely followed surface temperatures through 
0800 UTC on 22 January. Differences in surface temperatures 
between forecasts were relatively small and less than 1°C. 

Individual model forecasts exhibit larger differences when 
examining storm total precipitation (Fig. 3b), compared 
against observed precipitation that began around 0200 UTC 
on 22 January. The Goddard forecast was the first to pro- 
duce light precipitation at the CARE site, preceding ob- 
served accumulations by approximately four hours, while the 
WSM6, WDM6, Thompson, and SBU-Lin schemes lagged 
the observed precipitation onset by one to two hours. Since 
the double-moment version (WDM6) of the WSM6 retains 
the ice processes of the WSM6 and only provides a double 
moment representation for rain, accumulated precipitation in 
the two forecasts are equivalent, with no apparent impact 
from any upstream processes related to the rain category. 
All simulations follow the general trend in precipitation 
accumulation, but result in an underestimate of storm total 
accumulation through 0800 UTC when precipitation ended 
at the CARE site. All forecasts continued to accumulate 
precipitation beyond the observed ending time. The Morrison 
scheme, which includes double-moment representation of 
all precipitating species, obtained the minimum difference 
between simulated and accumulated precipitation ending 
at 0800 UTC and performed best overall when predicting 
hourly and storm total accumulations. 

6. Hydrometeor Profiles 

In a simulation of the 22 January 2007 event by Shi et al. 
(2010), the model forecast was deemed able to reproduce 
the general onset and character of precipitation. Molthan 
et al. (2010) demonstrated comparable precipitation coverage 
between the Goddard scheme forecast and radar observations 
at 0600 UTC, justifying comparisons between aircraft data 
and model profiles within 50 km of the King City radar. Here, 
the 50 km range of profiles is replicated, with conditional 
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Table 2. Characteristics of Microphysics Schemes Used in Generating WRF Model Forecasts 


Scheme 

Moments 

Notes 

Selected References 

Goddard 

1 

Saturation adjustment by Tao and Simpson (1993) 

Tao et al. (2003) 

WSM6 

1 

N os (T) by Houze et al. (1979) 

Hong et al. (2004) 

SBU-Lin 

1 

N os (T) by Houze et al. (1979) 

M(D) functions of diagnosed riming factor Ri,T 
V(D) functions of diagnosed riming factor Ri,T 

Lin and Colle (2010) 
Lin et al. (2010) 

Thompson 

1 

Predicts M n as fiMi-T) from Field et al. (2005) 
V(D) = a v D b 'e-f\f v = 125 

Thompson et al. (2008) 

WDM6 

2 

Double moment only applies to rain category 

Hong et al. (2010) 

Morrison 

2 

Number concentration and mass for each species 

Morrison et al. (2005) 


Table 3. Parameters Defining Relationships Within Microphysics Schemes Used in Generating WRF Model Forecasts 


Scheme 

Nos ( rn - 4 ) 

Ms 

p s (kg m~ 3 ) 

a m (kg m bm ) 

bm 

A^ (m 

a v ( m l bv s ! ) 

b v 

Goddard 

1.6x 10 7 

0 

100 

IPs 

3.0 

(6) 

1.305 

0.11 

WSM6 

f(T) 

0 

100 

IPs 

3.0 

(6) 

11.72 

0.41 

SBU-Lin 

f(T ) 

0 

m 

/( T,Rt) 

f(T,Ri) 

(6) 

f(T,Ri) 

f(T,Ri) 

Thompson 

N/A 

N/A 

m 

0.069 

2.0 

N/A 

40.0 

0.55 

WDM6 

f(T) 

0 

100 

IPs 

3.0 

(6) 

11.72 

0.41 

Morrison 

/(Mo, A,) 

0 

100 

IPs 

3.0 

(6) 

11.72 

0.41 


mean profiles of non-zero hydrometeor content shown in 
Fig. 4. Although these schemes produce up to six hydrom- 
eteor categories and their total ice contents are comparable 
to aircraft estimates available through CVI measurements, 
each partitions the total ice content in different ways based 
upon their varied assumptions and simulated processes. None 
of the forecasts produced an appreciable amount of graupel, 
in agreement with digital photographs of large, lightly rimed 
aggregates that occurred at the surface (Petersen et al. 2007). 
The Goddard scheme was the only forecast to produce an 
appreciable profile of cloud liquid water content, although 
no significant amounts of liquid water were detected by 
aircraft instrumentation. This is likely a result of the order 
of microphysics operations within the code and discussed 
within Molthan et al. (2010). 

Each of the forecasts differ substantially in their parti- 
tioning of total ice mass among the cloud ice and snow 
categories. In each forecast, cloud ice contributes to the 
eventual development of snow mass through autoconversion, 
accretion, and deposition. Differences in the handling of 
these processes lead to varied profiles of cloud ice and snow 
mass. Precipitation in the Goddard scheme is characterized 
by an upper level layer of cloud ice (4-6 km) which transi- 
tions to the snow category around 4 km. The SBU-Lin and 
Morrison forecasts also produced a layer of cloud ice from 4- 
6 km, but with a reduced mass content and a faster transition 
to the snow category. In the Thompson forecast, the vertical 
profile is dominated by the snow category. In the WSM6 
and WDM6 formulation, ice category profiles are equivalent, 
since double-moment representation applies only to the rain 
category. These schemes stand apart from other forecasts 
with a cloud ice profile that exceeds the snow category at 
altitudes of 1 km and above, representing precipitation as 


large number concentrations of pristine ice crystals rather 
than aggregates. 

7. Temperature and Water Vapor Profiles 

In addition to the solid or liquid species, each scheme 
handles the sources and sinks of water vapor through phase 
change processes, contributing to sources or sinks for the 
hydrometeor classes, and latent heating within the vertical 
profile. Mean temperature profiles were constructed for each 
scheme using the same sets of model vertical profiles ob- 
tained within 50 km of the King City radar installation. 
Absolute differences between the mean temperature profiles 
for each forecast were less than 0.5°, with the largest 
differences focused in the lowest 1-2 km of the vertical 
profile, and all forecasts exhibiting a slight warm bias in 
the entire vertical column (Fig. 5). 

In order to compare simulated water vapor profiles against 
aircraft data, simulated water vapor fields were converted to 
relative humidities with respect to water and ice, and reported 
as the maximum value at each model vertical level. Greater 
differences in relative humidity occur, given the variety of 
mechanisms for initiating precipitation within each scheme, 
and their respective methods for implementing saturation 
adjustments. Aircraft data indicate that the entire vertical 
column was saturated (supersaturated) with respect to water 
(ice), but each scheme obtains various levels of saturation 
depending upon their assumptions and parameterizations 
(Fig. 6). Molthan et al. (2010) noted the discrepancy in satu- 
ration levels that occurs within the Goddard scheme forecast 
at temperatures colder than -15°C, attributed to a temperature 
threshold within the Tao and Simpson (1993) saturation 
adjustment scheme, an assumption repeated here. The WSM6 
and WDM6 forecasts are not saturated (subsaturated) with 
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Fig. 2. Storm total, liquid equivalent precipitation through 0600 UTC on 22 January 2007 for each of simulations employed in this study. 




21/12 21/16 21/20 22/00 22/04 22/08 22/12 

force erst Hour 


Fig. 3. Comparisons of observed and predicted 2 m air temperatures 
at the CARE site (top) and storm total, liquid equivalent precipitation 
accumulation, based upon six WRF forecasts of the 22 January 2007 
snowfall event. 


respect to water (ice), related to their handling of depositional 
growth or saturation adjustment processes. The SBU-Lin 
forecast allows for a linear decrease in supersaturation with 
respect to ice between 500 and 300 hPa. Although the SBU- 
Lin forecast approaches saturation with respect to water in 


the lowest 1-2 km, values decrease with decreasing pressure. 
The Thompson and Morrison schemes produce the best 
representation of the water vapor profile, each maintaining 
values near saturation with respect to water except for a 
minor reduction near 4 km. 

8. Size Distribution Parameters 

Five of the six aforementioned schemes assign an expo- 
nential size distribution to precipitating hydrometeors, a spe- 
cial case of the gamma distribution (1) where the dispersion 
parameter ju s is set to zero. Aircraft estimates of N os and X s 
were acquired from C3VP aircraft profiles shown in Fig. 1, 
retaining parameters that provide a reliable fit (. R 2 >0.8) 
between the resulting exponential distribution and the actual 
PSD. Aircraft estimates of size distribution parameters can 
be compared against model assumptions or outputs for all 
forecasts except the Thompson scheme. Discussions related 
to the Thompson scheme are deferred to the next section 
where the performance of each scheme is examined in terms 
of distribution moments. 

Mean profiles of particle size distribution parameters were 
acquired from WRF model vertical profiles within 50 km of 
the King City radar, then compared against aircraft measure- 
ments (Fig. 7). The constant value of N os in the Goddard 
forecast was incapable of representing vertical variability 
in aircraft observations (Molthan et al. 2010). The WSM6, 
WDM6 and SBU-Lin schemes determine N os using Houze 
et al. (1979) function of temperature and provide for some 
variability in N os and X s with height, except for the lowest 
1-2 km of the vertical profile where temperatures are nearly 
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Fig. 4. Profiles of hydrometeor content obtained as conditional means from WRF model grid points within 50 km of the King City radar location, compared 
with observations of ice water content provided by the CVI instrument aboard the Convair-580 aircraft. Total ice profiles represent the sum of cloud ice, snow, 
and graupel categories, where available. Color coding of aircraft data represents the profiles described in Fig. 1. 



Temperature (°C) 


Fig. 5. Mean profiles of air temperature from WRF model grid points 
within 50 km of the King City radar location, compared with observations 
acquired fromthe Convair-580 aircraft data. 


isothermal. Values of a m and b m within the SBU-Lin scheme 
are determined from local calculation of the riming factor Ri 
(Lin and Colle 2010). These calculations of a m and b m are 
less than the | p s used within the WSM6 and WDM6 and 


permit a reduction in A* despite an increase in the simulated, 
snow mass content (Fig. 4). 

9. Moments of Particle Size Distributions 

The Thompson scheme is unique because it uses mi- 
crophysical processes to acquire ice mass (M 2 ) and then 
uses equations relating M 2 and temperature to calculate 
other moments needed in the simulation of microphysical 
processes. The remaining schemes use the exponential size 
distribution and parameters to determine remaining moments 
needed to simulate processes and rely upon either fixed 
values or temperature-dependent functions to characterize 
N os • Given the disparity in techniques used to characterize 
PSDs within the aforementioned schemes, comparisons are 
first examined here using profiles of various M n obtained as 
mean values from WRF model profiles within 50 km of the 
King City radar installation (Fig. 8). Analysis herein focuses 
on contributions from the snow mass category, since except 
the Morrison two-moment prediction assume a monodisperse 
assignment for the cloud ice class. Determinination of the 
size or number concentration of cloud ice crystals varies 
within each scheme but large number concentrations of very 
small particles were not found to contribute significantly to 
moments of order greater than Mq. 

Comparisons of Mo represent the ability to represent the 
total number concentration (Fig. 8a). Observations suggest 
that total number concentrations decrease between cloud top 
and cloud base as larger aggregates develop from mergers of 
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Fig. 6. Mean profiles of relative humidities with respect to water and ice from WRF model grid points within 50 km of the King City radar location, 
compared with observations acquired from the Convair-580 aircraft data. 




Fig. 7. Mean vertical profiles of N os and X s acquired from WRF model grid points within 50 km of the King City radar, based upon assumptions and 
predicted snow contents unique to each scheme. Model results are compared against PSDs and their parameter estimates, acquired from aircraft profiles shown 
in Fig. 1. 


smaller crystals, reducing the total number of particles. This 
coincides with increases in total ice water content through 
vapor diffusion and light riming. The Morrison scheme 
replicates the general trend in Mq decrease toward cloud 
base, likely benefitting from the double-moment represen- 
tation of both mass and total number concentration. The 
Thompson scheme forecast overestimated observed values 
of Mo and failed to represent the trend of increasing Mo 
with decreasing height, despite functional relationships of 
temperature and predicted M 2 that provide flexibility in size 
distribution assignment. Remaining single-moment schemes 


employ a variety of strategies for determining Mq. Since 
each scheme uses an exponential size distribution, Mo is the 
ratio between the distribution intercept and slope parameter. 
The SBU-Lin, WSM6 and WDM6 forecasts use the Houze 
et al. (1979) relationship for N os (T), and although they 
underestimate observations, they follow the general increase 
in observed N os with height (Fig. 7). As the lapse rate of 
temperature is reduced in the lowest 3 km (Fig. 5), variablity 
in similated N os decreases, while observed values continue to 
decrease due to aggregation. The ?l s parameter is calculated 
based upon the snow content and assignment of N os (6), 
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Fig. 8. Mean profiles of various moments of simulated particle size distributions for the snow category, based upon varying scheme assumptions, compared 
against moments derived from aircraft measured particle size distributions acquired in 5 s intervals. Color coding of aircraft data represents the profiles 
described in Fig. 1. 


with underestimates of N os contributing to underestimates 
of X s and particle size distributions with mean crystal sizes 
larger than suggested by aircraft PSDs. The Goddard forecast 
produced little variability in Mo within the lowest 4 km, 
attributable to the use of a constant size distribution intercept 
and snow bulk density. 

Determination of M2 is crucial for the Thompson scheme. 
Within the Thompson scheme, prediction of the total ice 
water content is related to M2 through the assignment of 
the mass-diameter relationship, and M2 is used to calcu- 
late additional moments, with M n as functions of M2 and 
temperature. This avoids the use of constants within the 
prescribed size distribution and allows for vertical variability 
observed in nature as a function of M2 and temperature (Field 
et al. 2005 ). In this particular event, the Thompson scheme 
underestimated M2, particularly in the lowest 3 km, where 
large aggregates were observed (Molthan et al. 2010 ). Due 
to the underestimate of M2 and the limited range in temper- 
ature within the observed and simulated profiles (Fig. 8c), 
predicted, higher order moments continued to underestimate 
aircraft observations. Attempts by the Thompson scheme 
to use M2 and temperature to predict other moments have 
merit, as mean profiles of Mi_ 6 replicate the general trend 
in vertical variabilty. However, it may be that the rapid onset 
of aggregation in this specific event is not well represented 
by the Field et al. ( 2005 ) relationships of f(M2,T) currently 
used within the scheme. Modification of the f(M2,T) rela- 
tionships to better fit this event may improve upon the current 


fit between aircraft and simulated PSD moments. 

Other moments shown in Fig. 8 demonstrate observed 
vertical trends and the ability of each scheme to replicate 
other characteristics of aircraft particle size distributions. 
Depending upon the assignment of PSD characteristics, 
successful representation of lower order moments do not 
guarantee successful estimate of higher order moments, since 
higher order moments are increasingly sensitive to weighting 
by relatively small number concentrations of large targets 
( 5 ). For example, although mean profiles of M\ acquired 
from the Morrison scheme underestimate aircraft observed 
values, this underestimation is reduced for higher order 
moments of M2 through M3, and the scheme produces mean 
profiles of M4 and M§ that provide a good representation 
of values within the ascending aircraft profile. In a single- 
moment representation, the SBU-Lin forecast slightly under- 
estimated M2 values throughout the bulk of the profile, but 
produced a reasonable depiction of M4 and M^, providing 
a better fit to aircraft observations than the two-moment 
representation of the Morrison scheme. Profiles of X s and 
N os for the SBU-Lin scheme indicate that both parameters 
are generally underestimated, but higher order moments are 
inversely proportional to X} +n ( 5 ), so that underestimates 
(overestimates) of X s (mean particle size) contribute to larger, 
predicted values of M^. Simulated decreases in X s and N os 
between cloud top and cloud base represent an ability to 
represent some of the effects of aggregation within the single- 
moment formulation, comparable to the Morrison double- 
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moment representation. The Thompson scheme is capable of 
representing the vertical trend in but would need some 
modification to internal functions of M 2 and temperature in 
order to increase predicted values of M^. For comparison, 
the WSM 6 and WDM 6 forecasts predict values of X s and N os 
that are too large, or a total particle number concentration 
and mean aggregate diameter less than observed, resulting in 
a substantial underestimate of M^. Mean profiles from the 
Goddard scheme struggle to represent the vertical variability 
observed in aircraft data due to strict adherence to constants, 
described previously by Molthan et al. (2010). 

10. Comparisons of Radar Reflectivity Profiles 

Radar reflectivity is often used as a postprocessed model 
variable since it quickly depicts the coverage, structure, 
and relative intensity of precipitation within the forecast 
domain. The methodology of Smith (1984) was used to 
convert predicted ice water content and size distribution to 
the simulated reflectivity value, based upon diameters of 
equivalent, pure ice spheres. The resulting equation for an 
equivalent radar reflectivity factor (Z e ) from forecast model 
output is: 


where the coefficient represents the ratio of the dielectric 
factors for water and ice, and a radar reflectivity factor Z 
based upon spheres with ice mass equivalent to simulated 
snow crystals. By including all possible variability in particle 
size distribution and mass-diameter relationships, Z can be 
obtained from forecast model output as: 

Z = (— ) f D 2b D^N(D)dD. ( 8 ) 

\ ft Pi J Zo 

The Goddard, WSM 6 , WDM 6 , and Morrison schemes 
represent snow crystals as spherical shapes within an expo- 
nential size distribution and with a fixed bulk density, so that 
a m = |p 5 , b m = 3, and /! = 0 (Table 3). This combination 
of variables results in Z proportional to and resulting 
reflectivity profiles will be comparable to the vertical profile 
of M$. The Thompson scheme uses a fixed a m and b m = 2, 
resulting in similarity to M4. Since the SBU-Lin forecast 
provides for flexibility in both a m and b m , moments used in 
the calculation of Z will vary between M 4 and M^. 

In order to evaluate model performance, the 0600 UTC 
volume scan of King City horizontally polarized radar re- 
flectivity is compared against values simulated from model 
output, based upon vertical profiles within 50 km of the radar 
location. Observed radar reflectivity was used to construct a 
contoured frequency with altitude diagram (CFAD, Yuter and 
Houze 1995) with increments of 2 dBZ and 500 m in order to 
demonstrate vertical variability (Fig. 9). The mean profile of 
simulated reflectivity provides a comparison between mod- 
eled values and the relative frequency of observations within 
the same altitude range. Differences in modeled and simu- 
lated reflectivity relate to the size distribution, mass-diameter 
relationship, and snow mass content for each scheme. 


Since the Thompson and SBU-Lin schemes use informa- 
tion from the ambient environment to allow for variabil- 
ity in PSD and density characteristics within the vertical, 
their mean profiles represent the general, vertical trend in 
modal reflectivity values shown within the King City CFAD. 
The WSM 6 /WDM 6 schemes incorporate a temperature- 
dependent particle size distribution and follow the vertical 
trend in radar reflectivity below 2 km, but the rapid transition 
from snow to cloud ice content contributes to a sharper 
decrease in radar reflectivity with altitude than observed by 
the King City radar. 

In this case, single-moment schemes incorporating vertical 
variability in snow characteristics are capable of representing 
the vertical trend in King City radar reflectivity, but the 
Thompson and SBU-Lin schemes provide a better fit by 
maintaining populations of large, precpitating ice crystals 
through a deeper portion of the vertical column. Although the 
Goddard scheme provides the best fit between CVI estimates 
of IWC and predicted snow mass, fixed values of p s and N os 
result in a reduced reflectivity lapse rate and a median profile 
that does not represent observed trends in the lowest 3 km. 
Above the 3 km level, the rapid transition between snow and 
cloud ice reduces the median reflectivity profile in a manner 
similar to the WSM 6 /WDM 6 forecasts. 
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Fig. 9. Mean vertical profiles of N os and X s acquired from WRF model grid 
points within 50 km of the King City radar, based upon assumptions and 
predicted snow contents unique to each scheme. Model results are compared 
against PSDs and their parameter estimates, acquired from aircraft profiles 
shown in Fig. 1. 


11. Comparisons of Terminal Fall Speed 
Relationships 

Terminal fall speeds of snow crystals were measured at 
the CARE site using a Hydrometeor Velocity and Shape 
Detector (HVSD, Barthazy et al. 2004), which images crys- 
tals passing between a series of digital detectors, and uses 
repeated imagers to estimate fall speeds. Resulting crystal fall 
speeds were provided by GyuWon Lee (McGill University) 
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after determining that wind did not bias fall speeds during 
the snowfall period. It is assumed herein that the reported 
fall velocities are the terminal velocities for each crystal. 
Observations of fall speed were binned for each HVSD 
maximum dimension bin diameter at 5 cm s -1 increments, 
accumulated for over 1 1 ,000 individual flakes observed from 
0200 to 0800 UTC (Fig. 10). Sizes and terminal velocities 
were fit to the power law form of Locatelli and Hobbs (1974), 
but were restricted to sizes greater than 1 mm and with at 
least 50 observations to account for limitations in HVSD 
detection of small particles and provide for a reasonable 
sample size. 

The best fit relationship for HVSD crystals produced 
values of a v = 110.083 cm l ~ bv s -1 and b v = 0.145, com- 
parable to Locatelli and Hobbs (1974) values for “unrimed 
radiating assemblages of dendries” ( a v = 80.0 cm l ~ bv s -1 , 
b v = 0.16). Molthan et al. (2010) attributed this similarity 
to the presence of large aggregates dominating the lowest 
levels of aircraft and dual-polarimetric radar observations. 
Simulated hydrometeor mass flux and resulting precipitation 
are sensitive to both the particle size distribution and the 
chosen relationship between diameter and fall speed. 

Comparisons between the HVSD best-fit relationship and 
each forecast scheme (Table 3) are shown in Fig. 10. The 
Goddard forecast underestimates fall speeds at all diameters, 
while the Morrison and WSM6/WDM6 forecasts overesti- 
mate (underestimate) fall speeds for particles larger than 
2 mm (smaller than 1 mm). The SBU-Lin scheme is the best 
fit to observations, although a slight overestimate occurs for 
all fall speeds, on the order of 0.1 m s -1 . The Thompson 
scheme oversetimates fall speeds for particles larger than 
1 mm, but includes an expoential decay term to reduce the 
fall speeds of large particles. Although the adjustment was 
not strong enough for this case, it is an improvement over 
the Morrison and WDM6AVSM6 where the overestimate 
continues to grow with increase in particle maximum dimen- 
sion. Although the Morrison scheme produced the best fit to 
observed surface precipitation accumulation (Fig. 3) and the 
simulation of aggregation benefits from the inclusion of the 
second moment (Fig. 7), the higher precipitation rate may 
have been obtained by overestimating particle fall speeds. 
Conversely, the SBU-Lin scheme provided the best overall 
fit to particle fall speeds, but prediced snow content was less 
than the mass acquired from CVI measurements (Fig. 4), 
and simulated particle size distributions may have produced a 
mean size larger than observations, combining for a precipita- 
tion rate closer to observations. In general, conclusions about 
model performance based upon sensible weather elements 
such as precipitation rate should also consider the reliability 
of assumptions in particle size distribution, mass-diameter 
relationship, and fall speeds. 

12. Summary and Conclusions 

An intensive observation period of the C3VP field cam- 
paign measured snow crystal particle size distributions, ice 
water content, and other atmospheric state variables within 
a broad region of snowfall associated with a passing midlat- 
itude cyclone on 22 January 2007. Simulations of the event 



Maximum Dimension (mm) 


Fig. 10. Mean vertical profiles of N os and X s acquired from WRF model 
grid points within 50 km of the King City radar, based upon assumptions and 
predicted snow contents unique to each scheme. Model results are compared 
against PSDs and their parameter estimates, acquired from aircraft profiles 
shown in Fig. 1. 


were performed with the Advanced Research WRF, using a 
baseline set of physical parameterizations with various selec- 
tions of single and double-moment microphysics schemes. 
The resulting model output fields of hydrometeor content, 
water vapor, and temperature were compared to in situ 
aircraft and surface measurements. Assumptions of particle 
size distribution, mass-diameter relationship, and diameter- 
terminal velocity relationships were evaluated using available 
C3VP datasets. 

Each of the single and double-moment schemes have vari- 
ous strengths and weaknesses, but all produced a reasonable 
simulation of the event, including surface temperatures and 
liquid equivalent precipitation rates. The representation of di- 
verse size distrbutions and particle effective bulk densities is 
best schieved by incorporating variability in size distribution 
parameters and mass-diameter relationships. These strategies 
are employed by the WSM6/WDM6, the Thompson, and 
the SBU-Lin forecasts which use a temperature-dependent 
size distribution intercept, temperature-dependent relation- 
ships between various size distribution moments, or the 
local prediction of size distribution, fall speed, and density 
characteristics as a function of crystal riming. Fixed values 
for these characteristics are often unable to represent changes 
in crystal characteristics throughout the vertical column. 
Two-moment representation assists in the depiction of the 
aggregation process by allowing for better representation of 
total particle number concentrations, but could be improved 
upon by allowing for greater flexibility in the remaining 
mass-diameter and fall speed relationships. It is unlikely 
and perhaps unrealistic to expect any given scheme to 
precisely simulate the characteristics of a single event, but 
field campaign data sets should be examined where available 
to evaluate the assumptions present within various physical 
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parameterizations given their increased utilization in current 
operational or experimental weather forecast models. In the 
interim, ensemble prediction strategies that combine various 
scheme outputs into a range of plausible events may assist 
in local, high resolution forecasts 
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transitioning unique NASA data and research technologies to operations 




Background 


■ Canadian CloudSat/CALIPSO Validation Project (C3VP) 

Designed to obtain surface and aircraft measurements of winter 
precipitation coincident with CloudSat/CALIPSO orbits. 

Participation by the NASA Global Precipitation Measurement 
mission ground validation component with emphasis on future 
development of precipitation retrievals. 

■ Field Campaign Observations 

■ Aircraft sampling of atmospheric state parameters, 
hydrometeor content, and particle size distributions. 

Standard surface measurements along with ground-based 
remote sensors, particle size distributions and fall speeds. 

■ Microphysics Schemes 

■ Prediction of one moment (mass) or two moments (mass and 
size distribution parameters) provided by a variety of previous 
studies, using C3VP data to evaluate performance. 



Snowfall Event 


■ 22 January 2007 

■ Widespread, moderate 
to heavy snowfall 
occurred throughout 
Southern Ontario 

■ Precipitation 
developed in advance 
of a mid-level trough 
and warm frontal 
boundary. 



Figure 1. Radar reflectivity from the C-band, 
dual-polarimetric radar at King City, Ontario, 
along with C3VP aircraft flight track and related 
observations. 



Model Configuration 


■ Use the WRF model in a 
configuration that can 
reproduce the onset and 
duration of the event. 

■ WRF V3.1, 24-hour forecast 
beginning at 12 UTC on 
January 21. 

■ Triple-nested grid with 
GFS 0.5 degree initial and 
boundary conditions 

■ Shi etal. 2010 
Molthan et al. 2010 

■ Model comparisons based 
upon output at 0600 UTC 
on 22 January 2007. 



Figure 2. Grid configuration for WRF model 
simulations of the 22 January 2007 snowfall 
event. 



Microphysics Schemes 


Overview of Snow Characteristics in Selected Microphysics Schemes 

Goddard 

(GSFC6G) 

• Single-moment, six class with graupel 

• Fixed distribution intercept 

• Fixed density spheres 

Tao et al. 2003 
Shietal. 2010 

WSM6 

• Single-moment 

• Size distribution intercept a function of temperature 

• Fixed density spheres 

Hong et al. 2005 

Thompson 

• Single-moment 

• Predicts snow content, then other moments are calculated and used 
within physics as a function of snow content and temperature. 

• Non-spherical mass-diameter relationship 

Thompson et al. 2008 

SBU-Lin 

• Single-moment 

• Uses a diagnosed riming factor, Ri, to characterize precipitating ice 

• Distribution intercept is a function of temperature 

• Predicts mass-diameter and diameter-fall speed characteristics from Ri 

Lin et al. 2010 
Lin and Colle 2010 

WDM6 

• Double moment rain category 

• Otherwise comparable to the WSM6 scheme for ice categories 

Hong et al. 2010 

Morrison 

• Double moment in all species 

• Fixed density spheres 

Morrison et al. 2005 




Hydrometeor Profiles 
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Figure 3. Conditional mean hydrometeor profiles acquired from model output within 50 km of the King City radar. 


Each scheme produces a unique separation of cloud ice and snow within the vertical 
column, and all tend to underestimate 1-3 km IWC, except for the Goddard scheme. 



WaterVapor Profiles 





Figure 4. Mean profiles of saturation with respect to ice and water for profiles within 50 km of the King City radar. 


Successful representation of the watervapor profile varies among schemes, based upon 
representation of depositional processes or saturation adjustment schemes. 



Size Distribution Parameters 
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Figure 5. Mean profiles of exponential size distribution parameters based upon scheme assumptions and predicted 
snow mass content, based upon profiles within 50 km of the King City radar. Aircraft estimates obtained from reliable 
(R 2 > 0.8) fits between exponential size distribution functions and measured particle size distributions. 


Fixed values of N os fail to represent natural variability, though N 0S (T) offers some 
improvement. Double-moment representations and non-spherical particle size 
distributions improve the fits for the distribution slope parameter, A s . 



Distribution Moments 
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Figure 6. As in Figure 5 but representing mean profiles of size distribution moments for the various schemes. 


M n = N os f(i+|j+n) / A (1+ ^ +n) or M n = f(M 2/ T) forThompson 


Double-moment representation improves profiles of number concentration (M 0 ). 
Temperature-dependent relationships of N os or M n provide a better representation of 
vertical variability, though they often underestimate observations, related to IWC. 



Terminal Velocity 


■ Terminal velocities 
obtained at the surface for 
observed aggregates. 

■ "HydrometeorVelocity and 
Shape Detector" for 
particles > 1 mm. 

Barthazyetal. 2004 

■ SBU-Lin scheme uses a 
riming factor that fits well 
to observations for this 
case. 

■ Other schemes over- or 
underestimate to varying 
degree. 
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Figure 7. Joint histogram of particle 
maximum dimension and fall speed, 
power-law fits from model forecasts, and 
C3VP relationship derived from mean fall 
speeds in each HVSD size bin. 
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Summary and Conclusions 


■ Observations from the C3VP storm of 22 January 2007 were used 
to evaluate various microphysics schemes within WRF V3.1. 

■ One event is insufficient to declare a "winner", but each scheme 
demonstrates some relative strengths and weaknesses: 

Fixed values for PSD parameters are insufficient. 

Temperature-based functions help, but could struggle in isothermal or 
inverted profiles. 

Vapor profiles are highly sensitive to assumed processes or saturation 
adjustment. 

Non-spherical, variable density particle assumptions and flexibility in 
characteristics often improves fits to observations. 

Fall speed relationships may impact IWC development and QPF. 

■ Ideally, relative strengths and weaknesses can be leveraged to 
make additional model improvements. 



Questions? 


■ Contact: andrew.molthan(a)nasa.qov 

■ SPoRT: http://weather.msfc.nasa.gov/sport 

■ Complete list of references is available within 
the extended abstract. 



